using System;
using System.Text;

/// <summary>
/// Utility functions for Matrices, and Probability Distributions.
/// </summary>
namespace Diversity.Maths
{
	/// <summary>
	/// Supports most standard matrix opperations on double based matrices
	/// </summary>
	public class Matrix
	{
	    private int _rows;
		private int _columns;
        private double[,] _matrix;

        /// <summary>
        /// The rows down the matrix
        /// </summary>
        public int Rows
        {
            get { return _rows; } 
        }

        /// <summary>
        /// The columns across the matrix
        /// </summary>
        public int Columns
        {
            get {return _columns;}
        }

        /// <summary>
        /// Elements of the matrix specified by row,column
        /// </summary>
        /// <param name="row">The number of rows</param>
        /// <param name="column">The number of columns</param>
        /// <returns></returns>
        public double this[int row, int column]
        {
            get { return _matrix[row, column];}
            set { _matrix[row, column] = value; }
        }

        /// <summary>
        /// Elements of the matrix specified by row,column
        /// </summary>
        /// <param name="rows">The number of rows</param>
        /// <param name="columns">The number of columns</param>
        public Matrix(int rows, int columns)
		{
			_rows = rows;
			_columns = columns;
            _matrix = new double[rows, columns];
		}

		#region Operators
 
        /// <summary>
        /// Equals.
        /// </summary>
        /// <param name="lmatrix">Left matrix</param>
        /// <param name="rmatrix">Right matrix</param>
        /// <returns>True if they are the same</returns>
		public static bool operator == (Matrix lmatrix, Matrix rmatrix)
		{
			return Equals(lmatrix, rmatrix);
		}

        public static bool operator !=(Matrix lmatrix, Matrix rmatrix)
		{
			return !(lmatrix == rmatrix);
		}

        public static Matrix operator *(Matrix leftMatrix, Matrix rightMatrix)
		{
			if(leftMatrix.Columns != rightMatrix.Rows)
				throw new Exception("leftMatrix columns must equal rightMatrix rows");

            Matrix result = new Matrix(leftMatrix.Rows, rightMatrix.Columns);

			for(int i = 0; i < leftMatrix.Rows; i++)
			{
				for(int j = 0; j < rightMatrix.Columns; j++)
				{
					for(int k = 0; k < rightMatrix.Rows; k++)
					{
						result[i,j] += leftMatrix[i,k] * rightMatrix[k,j];
					}
				}
			}

			return result;
		}

        public static Matrix operator *(Matrix matrix, double multiplier)
		{

            Matrix result = new Matrix(matrix.Rows, matrix.Columns);

			for(int i = 0; i < matrix.Rows; i++)
			{
				for(int j = 0; j < matrix.Columns; j++)
				{
					result[i,j] = matrix[i,j] * multiplier;
				}
			}

			return result;
		}

        public static Matrix operator *(double multiplier, Matrix matrix)
		{
			return matrix * multiplier;
		}

        public static Matrix operator /(Matrix matrix, double divisor)
		{
            return matrix*(1.0/divisor);
		}

        public static Matrix operator +(Matrix leftMatrix, Matrix rightMatrix)
		{

			if(leftMatrix.Rows != rightMatrix.Rows || leftMatrix.Columns != rightMatrix.Columns)
				throw new Exception("Plus operator requires matrices of the same size");

            Matrix result = new Matrix(leftMatrix.Rows, leftMatrix.Columns);

			for(int i = 0; i < leftMatrix.Rows; i++)
			{
				for(int j = 0; j < leftMatrix.Columns; j++)
				{
					result[i,j] = leftMatrix[i,j] + rightMatrix[i,j];
				}
			}

			return result;
		}

        public static Matrix operator -(Matrix leftMatrix, Matrix rightMatrix)
		{

			if(leftMatrix.Rows != rightMatrix.Rows || leftMatrix.Columns != rightMatrix.Columns)
                throw new Exception("Minus operator requires matrices of the same size");

            Matrix result = new Matrix(leftMatrix.Rows, leftMatrix.Columns);

			for(int i = 0; i < leftMatrix.Rows; i++)
			{
				for(int j = 0; j < leftMatrix.Columns; j++)
				{
					result[i,j] = leftMatrix[i,j] - rightMatrix[i,j];
				}
			}

			return result;
		}

		/// <summary>
		/// Returns the transpose of the matrix
		/// i.e. rows become columns
		/// </summary>
		/// <param name="matrix"></param>
		/// <returns>The transpose of the matrix</returns>
        public static Matrix operator ~(Matrix matrix)
		{

            Matrix result = new Matrix(matrix.Columns, matrix.Rows);

			for(int i = 0; i < matrix.Rows; i++)
			{
				for(int j = 0; j < matrix.Columns; j++)
				{
					result[j,i] = matrix[i,j];
				}
			}

			return result;
		}


        public static Matrix operator !(Matrix matrix)
		{
			if(matrix.Determinant() == 0)
				throw new Exception("Cannot invert a singular matrix");

			if(matrix.Rows == 2 && matrix.Columns == 2)
			{
                Matrix tempMatrix = new Matrix(2, 2);

				tempMatrix[0,0] = matrix[1,1];
				tempMatrix[0,1] = -matrix[0,1];
				tempMatrix[1,0] = -matrix[1,0];
				tempMatrix[1,1] = matrix[0,0];

				return tempMatrix / matrix.Determinant();
			}
			return matrix.Adjugate()/matrix.Determinant();

		}

		#endregion

		#region Methods



		public double Determinant()
		{
            double determinent = 0;

			if(_rows != _columns)
				throw new Exception("Matrix isn't square");

            if (_rows == 2 && _columns == 2)
            {
                determinent = (_matrix[0, 0] * _matrix[1, 1]) - (_matrix[0, 1] * _matrix[1, 0]);
                return determinent;
            }
            else
            {
                for (int j = 0; j < _columns; j++)
                {
                    Matrix tempMatrix = Minor(0, j);
                    determinent += (int) Math.Pow(-1, j)*this[0, j]*tempMatrix.Determinant();
                }
                return determinent;
            }
		}

        public Matrix Adjugate()
		{
			if(_rows < 2 || _columns < 2)
                throw new Exception("Only square matrixes have adjugates");

            Matrix adjugateMatrix = new Matrix(_columns, _rows);

			for(int i = 0; i < _rows; i++)
			{
				for(int j = 0; j < _columns;j++)
				{
                    Matrix tempMatrix = Minor(i, j);
					adjugateMatrix[j,i] = (int)Math.Pow(-1,i+j) * tempMatrix.Determinant();
				}
			}
			return adjugateMatrix;
		}

        public Matrix Minor(int row, int column)
		{
			if(_rows < 2 || _columns < 2)
				throw new Exception("Minor not available");

            Matrix minorMatrix = new Matrix(_rows - 1, _columns - 1);
 
			for(int k = 0; k < minorMatrix.Rows; k++)
			{
                int i;
				if(k >= row)
					i = k + 1;
				else
					i = k;

				for(int l = 0; l < minorMatrix.Columns; l++)
				{
                    int j;
					if(l >= column)
						j = l + 1;
					else
						j = l;

					minorMatrix[k,l] = this[i,j];
				}
			}

			return minorMatrix;
		}

		public bool IsInvertible()
		{
			if(Determinant() == 0)
			{
				return false;
			}
			return true;
		}

        public static Matrix Identity(int size)
        {
            Matrix identityMatrix = new Matrix(size,size);
            for (int i=0; i < size; i++)
            {
                identityMatrix[i, i] = 1;
            }
            return identityMatrix;
        }

        #endregion

        #region Overrides

        public override string ToString()
        {
            StringBuilder sb = new StringBuilder();

            for (int i = 0; i < _rows; i++)
            {
                for (int j = 0; j < _columns; j++)
                {
                    sb.Append(_matrix[i, j].ToString("0.0000")).Append("\t\t");
                }
                sb.Append("\n");
            }
            return sb.ToString();
        }

        public override bool Equals(Object obj)
        {
            Matrix matrix = (Matrix)obj;

            if (_rows != matrix.Rows || _columns != matrix.Columns)
                return false;

            for (int i = 0; i < _rows; i++)
            {
                for (int j = 0; j < _columns; j++)
                {
                    if (this[i, j] != matrix[i, j])
                        return false;
                }
            }

            return true;
        }

        public override int GetHashCode()
        {
            double result = 0;

            for (int i = 0; i < _rows; i++)
            {
                for (int j = 0; j < _columns; j++)
                {
                    result += this[i, j];
                }
            }
            return (int)result;
        }

       
        #endregion

	}
}
